Introducción a R como GIS

## Loading required package: sp
## Checking rgeos availability: TRUE

Resumen

En este tutorial vermos como utilizar R como un Sistema de Información Geográfica (SIG). R utiliza un conjunto de librerías para procesar objetos espaciales; las liberías que vamos a utilizar en este tutorial son:

  • sp: Classes and methods for spatial data; the classes document where the spatial location information resides, for 2D or 3D data. Utility functions are provided, e.g. for plotting data as maps, spatial selection, as well as methods for retrieving coordinates, for subsetting, print, summary, etc.
  • rgdal: Provides bindings to Frank Warmerdam’s Geospatial Data Abstraction Library (GDAL) (>= 1.6.3) and access to projection/transformation operations from the PROJ.4 library. The GDAL and PROJ.4 libraries are external to the package, and, when installing the package from source, must be correctly installed first. Both GDAL raster and OGR vector map data can be imported into R, and GDAL raster data and OGR vector data exported. Use is made of classes defined in the sp package. Windows and Mac Intel OS X binaries (including GDAL, PROJ.4 and Expat) are provided on CRAN.
  • raster: Reading, writing, manipulating, analyzing and modeling of gridded spatial data. The package implements basic and high-level functions. Processing of very large files is supported.
  • rgeos: Interface to Geometry Engine - Open Source (GEOS) using the C API for topology operations on geometries. The GEOS library is external to the package, and, when installing the package from source, must be correctly installed first. Windows and Mac Intel OS X binaries are provided on CRAN.
  • leaflet: Create and customize interactive maps using the ‘Leaflet’ JavaScript library and the ‘htmlwidgets’ package. These maps can be used directly from the R console, from ‘RStudio’, in Shiny apps and R Markdown documents.
  • rasterVis: Methods for enhanced visualization and interaction with raster data. It implements visualization methods for quantitative data and categorical data, both for univariate and multivariate rasters. It also provides methods to display spatiotemporal rasters, and vector fields. See the website for examples.

En general trabajaremos sobre tareas específicas y vermos como ejecutarlas desde R. Finalmente se asume que los estudiantes ya tiene conocimientos básicos de SIG como:

Vector
  • Qué es un vector
  • Cuantos tipos de vectores hay?
    • Puntos
    • Líneas
    • Polígonos
  • Conocer la estructura interna de un vector
    • Número de vectores
    • Tabla de atributos
Raster
  • qué es un raster?
  • Qué es un extent?
  • Cuáles son sus dimensiones?
  • Saber qué es una proyección?
  • Saber si un raster está proyectado o no?

Asignar el directorio de trabajo

Trabajando con vectores desde R

Como se mencionó existen tre tipos de vectores espaciales; estos son puntos, líneas y polígonos. R utiliza un conjunto de librerías para procesar objetos espaciales. Las librerías que utilizamos para leer y hacer operaciones sobre vectores son: rgdal y sp

Cómo leo un vector desde R.

El comando para leer un shapefile es readOGR; este se encuentra en la librería rgdal. Los argumentos del comando son:

  • dsn: Es el folder donde tenemos nuestro archivo vector (shapefile)
  • layer: El nombre de la capa que queremos leer
## rgdal: version: 1.3-6, (SVN revision 773)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.1.3, released 2017/20/01
##  Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rgdal/gdal
##  GDAL binary built with GEOS: FALSE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rgdal/proj
##  Linking to sp version: 1.3-1
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/Luis/Dropbox/CursoNichos2019/destdv1gw", layer: "destdv1gw"
## with 382 features
## It has 8 fields
## Integer64 fields read as strings:  COV_ COV_ID

Grafiquemos el mapa

Explorando los atributos del vector

Uno de los comandos de más útilies para explorar los elementos que conforman a un objerto, es el comando str.

o bien para darnos solo una idea de su estrutura podemos escribir el nombre del objeto

Los objetos espaciales en R son de la clase S4; para extraer los atributos de estos objetos (indexar) se utiliza el símbolo de **@**. Veamos como extraer la tabla de atributos del mapa map_vec

AREA PERIMETER COV_ COV_ID ENTIDAD CAPITAL RASGO_GEOG CVE_EDO
0 6.7179224 18.9211577 2 1 BAJA CALIFORNIA Mexicali NA 02
1 16.6952450 30.6009734 3 5 SONORA Hermosillo NA 26
2 0.0000470 0.0318774 4 2 BAJA CALIFORNIA NA ISLA 02
3 0.0000177 0.0175657 5 3 BAJA CALIFORNIA NA ISLA 02
4 0.0001608 0.0679551 6 4 BAJA CALIFORNIA NA ISLA 02
5 0.0127180 0.5160367 7 8 BAJA CALIFORNIA NA ISLA 02

Número de vectores o poligonos

## [1] 382

Guardar un shapefile

Ahora guardamos como un shapefile separado, el polígono

Filtrar los puntos que caen dentro de un polígono

Generemos unos puntos aleatorios definidos en el extent del mapa; el comando nombre_objeto@bbox nos muestra el extent de nuestro shapefile.

##          min       max
## x -118.36603 -86.71074
## y   14.53401  32.71877

Extraer los atributos en los puntos del poligono

longitude latitude AREA PERIMETER COV_ COV_ID ENTIDAD CAPITAL RASGO_GEOG CVE_EDO
2 -110.09291 29.04229 16.695245 30.60097 3 5 SONORA Hermosillo NA 26
10 -97.82715 19.84400 2.939128 18.83949 298 313 PUEBLA Herorica Puebla de Zaragoza NA 21
14 -97.09354 19.36795 6.062359 33.15909 255 269 VERACRUZ DE IGNACIO DE LA LLAVE Xalapa-Enriquez NA 30
20 -103.97467 27.88340 22.882945 29.38276 10 10 CHIHUAHUA Chihuahua NA 08
22 -98.87711 23.35904 6.917870 30.18376 71 75 TAMAULIPAS Ciudad Victoria NA 28
24 -93.17818 16.89131 6.214071 17.55246 371 389 CHIAPAS Tuxtla Gutierrez NA 07

Reproyectar un polígono

Con R podemos hacer reproyección de coordenadas, podemos pasar de coordenadas geográficas (logitude, latitud) a coordenadas utm (metrosal nivel del mar).

Veamo como proyectar nuesto mapa de Puebla a coordenadas utm

Regresemos a coordenadas geográficas

Descargar datos espaciales desde R

Ahora veremos como descargar datos espaciales desde R; para ello utlizaremos la función getData del paquete raster. La función permite descargar datos de diferentes bases de datos como:

Límites territoriales

  • GADM: Datos de límites territoriales de países

Datos climáticos

Datos Digitales de Elevación

  • SRTM: Datos Digitales de Elevación SRTM 90m.
  • mapzen: Datos Digitales de Elevación aprox ~ 9m.

Trabajando con rasters desde R

La librería estrella para trabjar con capas raster es raster. Esta es un wrapper que se comunica con funciones de alto nivel escritas en el lenguaje C para procesar objetos raster. Es importante notar que los raster no son nada más que arreglos bidimensionales de tipo numérico o categórico.

Al igual que los archivos de tipo vector tienen un extent, sin embargo tienen otros atributos como dimensión y resolución espacial. Entre mayor sea la resolución del raster, mayor será su dimensión (mayor número de filas y columnas).

Cómo leer un raster

El comando para leer 1 solo raster es raster

Cómo leer un conjunto (stack) de capas

A veces (casi siempre) será necesario leer un conjunto de capas que tienen la misma resolución y extent (i.e. WorldClim layers). El comando stack de raster permite crear stack de capas; para ello es necesario indicar con en un vector los path o rutas a la carpeta que contiene nuestras capas.

Veamos como leer o almacenar en la memoria de R el conjunto de capas de WorldClim layers

Primero en guardamos en una variable los path a nustras capas

Usamos el comando stack

Graficamos las primer 5 capas

notemos que utilice la misma sintaxis que utilizamos para indexa listas…

Hacer un crop

La función crop utiliza el extent de un shapfile para hacer recortar un raster o un stack. Veamos un recorte de las capas de WorldClim con el poligono de puebla.

Hacer una mascara

Al igual que crop, mask hace un recorte pero conserva la forma del poligono.

Guardar un raster

El comando para gurardar un raster es wirteRaster y está en el paquete raster

Convertir a puntos un stack

x y bio1 bio10 bio11 bio12 bio13 bio14 bio15 bio16 bio17 bio18 bio19 bio2 bio3 bio4 bio5 bio6 bio7 bio8 bio9
-97.75000 20.75000 240 274 192 1265 202 37 56 561 125 522 133 106 53 3329 332 134 198 270 202
-97.75000 20.58333 237 270 189 1737 282 49 64 811 150 784 165 109 54 3287 330 130 200 266 199
-97.91667 20.41667 222 254 178 2278 421 48 77 1179 151 905 161 119 56 3054 321 111 210 248 188
-98.25000 20.25000 148 171 122 1219 246 23 77 630 81 298 115 134 65 1923 249 45 204 157 133
-98.08333 20.25000 176 203 141 2142 402 42 79 1144 144 480 149 128 60 2506 277 67 210 194 151
-97.91667 20.25000 199 228 159 2587 480 55 77 1359 180 985 190 121 58 2794 298 91 207 221 169

Extraer valores de un raster o un stack

Podemos extraer los valores de un stack utilizando coordenadas.

bio1 bio10 bio11 bio12 bio13 bio14 bio15 bio16 bio17 bio18 bio19 bio2 bio3 bio4 bio5 bio6 bio7 bio8 bio9
2 231 293 161 465 144 3 110 313 19 279 80 158 51 5217 381 75 306 288 225
9 167 256 74 271 54 6 71 135 20 122 26 160 45 7076 338 -13 351 245 125
10 171 196 137 1832 372 42 70 884 144 413 158 124 61 2358 267 66 201 187 147
13 221 290 139 639 119 23 50 242 92 182 98 130 43 5880 362 64 298 265 150
14 149 169 124 1721 321 41 76 872 127 460 160 120 66 1800 238 57 181 160 133
15 175 259 82 683 88 28 37 235 96 212 96 145 42 6877 333 -5 338 225 82

Tabla de correlaciones…

bio1 bio10 bio11 bio12 bio13 bio14 bio15 bio16 bio17 bio18 bio19 bio2 bio3 bio4 bio5 bio6 bio7 bio8 bio9
bio1 1.00 0.65 0.87 0.35 0.41 0.07 0.08 0.39 0.07 0.33 0.14 -0.52 0.28 -0.44 0.45 0.88 -0.60 0.56 0.68
bio10 0.65 1.00 0.20 -0.02 -0.14 0.21 -0.33 -0.15 0.21 -0.10 0.19 -0.34 -0.48 0.40 0.88 0.26 0.16 0.49 0.34
bio11 0.87 0.20 1.00 0.44 0.60 -0.06 0.33 0.60 -0.07 0.47 0.04 -0.43 0.68 -0.82 0.03 0.98 -0.87 0.40 0.69
bio12 0.35 -0.02 0.44 1.00 0.90 0.71 -0.32 0.91 0.70 0.82 0.77 -0.66 0.22 -0.43 -0.30 0.55 -0.62 -0.09 0.40
bio13 0.41 -0.14 0.60 0.90 1.00 0.37 0.08 0.99 0.36 0.83 0.45 -0.50 0.49 -0.65 -0.35 0.64 -0.73 0.04 0.43
bio14 0.07 0.21 -0.06 0.71 0.37 1.00 -0.81 0.38 0.99 0.51 0.94 -0.64 -0.39 0.17 -0.08 0.10 -0.12 -0.27 0.18
bio15 0.08 -0.33 0.33 -0.32 0.08 -0.81 1.00 0.08 -0.82 -0.06 -0.70 0.51 0.67 -0.50 -0.11 0.17 -0.20 0.25 0.07
bio16 0.39 -0.15 0.60 0.91 0.99 0.38 0.08 1.00 0.37 0.82 0.47 -0.48 0.49 -0.65 -0.35 0.63 -0.72 0.02 0.43
bio17 0.07 0.21 -0.07 0.70 0.36 0.99 -0.82 0.37 1.00 0.51 0.95 -0.64 -0.41 0.18 -0.08 0.09 -0.11 -0.25 0.18
bio18 0.33 -0.10 0.47 0.82 0.83 0.51 -0.06 0.82 0.51 1.00 0.59 -0.53 0.32 -0.49 -0.37 0.52 -0.63 0.03 0.36
bio19 0.14 0.19 0.04 0.77 0.45 0.94 -0.70 0.47 0.95 0.59 1.00 -0.62 -0.28 0.07 -0.09 0.19 -0.21 -0.27 0.33
bio2 -0.52 -0.34 -0.43 -0.66 -0.50 -0.64 0.51 -0.48 -0.64 -0.53 -0.62 1.00 0.13 0.21 0.12 -0.61 0.60 -0.16 -0.36
bio3 0.28 -0.48 0.68 0.22 0.49 -0.39 0.67 0.49 -0.41 0.32 -0.28 0.13 1.00 -0.92 -0.42 0.57 -0.70 0.03 0.35
bio4 -0.44 0.40 -0.82 -0.43 -0.65 0.17 -0.50 -0.65 0.18 -0.49 0.07 0.21 -0.92 1.00 0.49 -0.77 0.90 -0.08 -0.45
bio5 0.45 0.88 0.03 -0.30 -0.35 -0.08 -0.11 -0.35 -0.08 -0.37 -0.09 0.12 -0.42 0.49 1.00 0.00 0.44 0.43 0.20
bio6 0.88 0.26 0.98 0.55 0.64 0.10 0.17 0.63 0.09 0.52 0.19 -0.61 0.57 -0.77 0.00 1.00 -0.90 0.38 0.71
bio7 -0.60 0.16 -0.87 -0.62 -0.73 -0.12 -0.20 -0.72 -0.11 -0.63 -0.21 0.60 -0.70 0.90 0.44 -0.90 1.00 -0.15 -0.55
bio8 0.56 0.49 0.40 -0.09 0.04 -0.27 0.25 0.02 -0.25 0.03 -0.27 -0.16 0.03 -0.08 0.43 0.38 -0.15 1.00 0.09
bio9 0.68 0.34 0.69 0.40 0.43 0.18 0.07 0.43 0.18 0.36 0.33 -0.36 0.35 -0.45 0.20 0.71 -0.55 0.09 1.00

Reproyectar un raster

Al igual que con los vectores, podemos reproyectar a los rasters. Lo anterior se puede hacer mediante el comando projectRaster de raster

Regresemos a geográficas

La belleza de lo interactivo

En este ejemplo veremos como limpiar datos de precias de una especie con un mapa interactivo

Leemos los datos

name longitude latitude prov date key
Puma concolor -80.68902 25.57658 gbif 2019-01-02 1986571407
Puma concolor -108.56825 33.12471 gbif 2019-01-04 1986581111
Puma concolor -116.57024 32.74562 gbif 2019-01-05 1986595849
Puma concolor -93.47792 16.05619 gbif 2018-01-08 1831091833
Puma concolor -120.73842 35.32380 gbif 2018-01-30 1837077255
Puma concolor -120.25755 39.43723 gbif 2018-01-24 1807289765
  1. Quitamos datos duplicados

Visualizamos los datos

Limpiamos datos

Escribimos nuestra tabla de datos

Descarga de datos de presencia spocc

Un poco más profesional

## class       : Extent 
## xmin        : -118.4042 
## xmax        : -86.7014 
## ymin        : 14.55055 
## ymax        : 32.71846

Aún más pro

Un ejemplo de programación funcional con purrr

  • Leer la base de datos de reptiles y anfibios
    • Partir el data.frame por especie.
    • Crear polígnos mínimos convexos de los registros.
    • Cortar las capas de WorldClim con el polígono.
    • Escribir los rasters (capas de modelación).
    • Buscar los datos de gbif usando polígonos.
    • Limpiar los duplicados espaciales.
    • Regresar el raster y los puntos de gbif.
## 
## Attaching package: 'magrittr'
## The following object is masked from 'package:raster':
## 
##     extract

Definimos la función para crear y escribir polígonos convexos en el ambiente global

Crando la función

Veamos con una especie

Grafiquemos

Ahora hagámoslo con algunas (ustedes pueden con todas las sps)

Exploremos la lista

Revisen su directorio de trabajo

Luis Osorio-Olvera ()